\documentclass[a4paper]{article}
\usepackage{amsmath,amssymb,amsfonts,latexsym,graphicx}
\usepackage[boldfont,slantfont]{xeCJK}
\usepackage{xcolor}
\usepackage{enumerate}
\setCJKmainfont{AR PL UMing CN}


\parindent=0mm
\textheight=20cm \textwidth=14cm
\parskip=0mm
\renewcommand\baselinestretch{1.2}
\oddsidemargin30pt \evensidemargin30pt


\begin{document}

\baselineskip 18pt
\parskip 10pt
\parindent 0em

\def\chntoday{\the\year~年~\the\month~月~\the\day~日}

\renewcommand\figurename{图}\renewcommand\tablename{表}
\renewcommand{\refname}{参考文献}
\renewcommand\abstractname{\bf\large 摘~~要}

\title{计算流体力学课程上机作业(一)}
\author{叶时炜}
\date{\chntoday}
\maketitle

\begin{abstract}
对线性气动方程，用迎风格式和Lax-Wendroff格式求不同初值条件下的程序。
\end{abstract}

\section{问题}

已知气动方程：
\begin{displaymath}\label{eq:qidong}
  \mathbf{U}_t+\mathcal{A} \mathbf{U}=0
\end{displaymath}
其中
\begin{displaymath}
  \mathbf{U}=(\rho,u)^T
\end{displaymath}

\begin{displaymath}
  \mathcal{A}=\left( 
  \begin{array}{cc}
    0&\rho_0\\
    \frac{a^2}{\rho_0}&0\\
  \end{array}\right)
\end{displaymath}

$\rho_0,a$均为常数。
\paragraph \space 编制迎风格式，Lax-wendroff格式的程序。要求初值$\mathbf{U}(x,0)$取两种情形：

\begin{enumerate}[(a)]
\item 光滑的周期函数。
\item $\mathbf{U}(x,0)$取为Riemann数据。
\end{enumerate}

\section{数值方法}
解出$\mathbf{A}$的特征值$\lambda_1=a,\lambda_2=-a$对应特征向量分别为：
\begin{displaymath}
  \begin{array}{cc}
    R(1)=\left( \begin{array}{c}
      \rho_0\\
      a\\
    \end{array}\right)&
    R(2)=\left(\begin{array}{c}
      -\rho_0\\
      a\\
    \end{array}\right)\\
  \end{array}
\end{displaymath}
设：
\begin{displaymath}
  \begin{array}{cc}
    R=\left( \begin{array}{cc}
      \rho_0&-\rho_0\\
      a&a\\
    \end{array}\right)&
    \mathbf{W}=R^{-1}\mathbf{U}\\
  \end{array}
\end{displaymath}
于是原方程变为：
\begin{displaymath}\label{eq:qidong}
  \mathbf{W}_t+\mathcal{B} \mathbf{W}_x=0
\end{displaymath}
其中
\begin{displaymath}
\mathcal{B}=R^{-1}\mathcal{A}R=\left( 
  \begin{array}{cc}
    a&0\\
    0&-a\\
  \end{array}\right)
\end{displaymath}
于是：
\begin{displaymath}
\left\{
  \begin{array}{c}
    \mathbf{W}(1)_t+a\mathbf{W}(1)_x=0\\
    \mathbf{W}(2)_t-a\mathbf{W}(2)_x=0\\
  \end{array}
\right.
\end{displaymath}
然后用上式，差分求解W（x，t）。利用$\mathbf{U}=R\mathbf{W}$得到$\mathbf{U}$

\section{数值实验}

简单起见，设$a=1$,$\rho_0=1$,相应的有
\begin{displaymath}
\begin{array}{cc}
    R=\left( \begin{array}{cc}
      1&-1\\
      1&1\\
    \end{array}\right) &
    R^{-1}=\left( \begin{array}{cc}
      \frac{1}{2}&\frac{1}{2}\\
      -\frac{1}{2}&\frac{1}{2}\\
    \end{array}\right)\\
\end{array}
\end{displaymath}

\subsection{周期边界条件}
先考虑周期边界:
\begin{displaymath}
\left\{\begin{array}{c}
\rho(x,0)=2+sin(\pi x)\\
u(x,0)=2+cos(\pi x)\\
\end{array}\right.
\end{displaymath}
于是有：
\begin{displaymath}
\left\{\begin{array}{c}
W1(x,0)=\frac{1}{2}(4+cos(\pi x)+sin(\pi x))\\
W2(x,0)=\frac{1}{2}(cos(\pi x)-sin(\pi x))\\
\end{array}\right.
\end{displaymath}
易得其精确解为：
\begin{displaymath}
\left\{\begin{array}{c}
W1(x,t)=\frac{1}{2}(4+cos(\pi(x-t))+sin(\pi(x-t)))\\
W2(x,t)=\frac{1}{2}(cos(\pi(x+t))-sin(\pi(x+t)))\\
\end{array}\right.
\end{displaymath}
于是相对应的$U$的精确解为：
\begin{displaymath}
\left\{\begin{array}{c}
\rho(x,t)=2+cos(\pi t)sin(\pi x)+sin(\pi t)sin(\pi x)\\
u(x,t)=2+cos(\pi t)cos(\pi x)-sin(\pi t)cos(\pi x)\\
\end{array}\right.
\end{displaymath}

其中$\rho$的等高线图像类似图\ref{ppr}
\begin{figure}[h!]\centering
\includegraphics[width=14cm]{../Periodic-Precise-rou.eps}
\caption{周期初始条件的精确解$\rho$}\label{ppr}
\end{figure}


\subsection{Riemann边界条件}
考虑Riemann数据:
\begin{displaymath}
\rho(x,0)=\left\{\begin{array}{cc}
1&x\leq0\\
3&x>0\\
\end{array}\right.
\end{displaymath}
\begin{displaymath}
u(x,0)=\left\{\begin{array}{cc}
5&x\leq0\\
1&x>0\\
\end{array}\right.
\end{displaymath}

于是有：
\begin{displaymath}
W1(x,0)=\left\{\begin{array}{cc}
3&x\leq0\\
2&x>0\\
\end{array}\right.
\end{displaymath}
\begin{displaymath}
W2(x,0)=\left\{\begin{array}{cc}
2&x\leq0\\
-1&x>0\\
\end{array}\right.
\end{displaymath}
易得其精确解为：

\begin{displaymath}
W1(x,t)=\left\{\begin{array}{cc}
3&x\leq t\\
2&x>t\\
\end{array}\right.
\end{displaymath}
\begin{displaymath}
W2(x,0)=\left\{\begin{array}{cc}
2&x\leq -t\\
-1&x>-t\\
\end{array}\right.
\end{displaymath}
于是相对应的$U$的精确解为：

\begin{displaymath}
\rho(x,t)=\left\{\begin{array}{cc}
1&x\leq -t\\
4&-t\leq x <t\\
3&x>t\\
\end{array}\right.
\end{displaymath}
\begin{displaymath}
W2(x,0)=\left\{\begin{array}{cc}
5&x\leq -t\\
2&-t \leq x < t\\
1&x\geq t\\
\end{array}\right.
\end{displaymath}

其中$\rho$的等高线图像类似图\ref{rpr}
\begin{figure}[h!]\centering
\includegraphics[width=14cm]{../Riemann-Precise-rou.eps}
\caption{Riemann初始条件的精确解$\rho$}\label{rpr}
\end{figure}

\section{实验结果}

\subsection{迎风格式}

\begin{enumerate}
\item 周期初始条件，$\rho$的计算结果等高线图像类似图\ref{pur}
\begin{figure}[h!]\centering
\includegraphics[width=14cm]{../Periodic-upwind-rou.eps}
\caption{周期初始条件的迎风格式计算解$\rho$}\label{pur}
\end{figure}
\item Riemann初始条件，$\rho$的计算结果等高线图像类似图\ref{rur}
\begin{figure}[h!]\centering
\includegraphics[width=14cm]{../Riemann-upwind-rou.eps}
\caption{Riemann初始条件的迎风格式计算解$\rho$}\label{rur}
\end{figure}
\end{enumerate}

收敛阶分析，对x轴\[-1,1\]分为40，80，160，320个网格计算了误差，作出下图。
\begin{enumerate}
\item $L^1$收敛阶，如图\ref{l1u}，双对数坐标下。
\begin{figure}[h!]\centering
\includegraphics[width=14cm]{../L1-convegence-order-of-upwind.eps}
\caption{网格大小与$L^1$误差的双对数图}\label{l1u}
\end{figure}
\item $L^2$收敛阶，如图\ref{l2u},双对数坐标。
\begin{figure}[h!]\centering
\includegraphics[width=14cm]{../L2-convegence-order-of-upwind.eps}
\caption{网格大小与$L^2$误差的双对数图}\label{l2u}
\end{figure}
\item $L^{\infty}$收敛阶，如图\ref{liu},双对数坐标。
\begin{figure}[h!]\centering
\includegraphics[width=14cm]{../L-infinity-convegence-order-of-upwind.eps}
\caption{网格大小与$L^{\infty}$误差的双对数图}\label{liu}
\end{figure}

\end{enumerate}
由图可知，迎风格式是一阶收敛。

\subsection{Lax-Wendrof格式}

\begin{enumerate}
\item 周期初始条件，$\rho$的计算结果等高线图像类似图\ref{pur}
\begin{figure}[h!]\centering
\includegraphics[width=14cm]{../Periodic-Laxwendrof-rou.eps}
\caption{周期初始条件的Lax-Wendrof格式计算解$\rho$}\label{pur}
\end{figure}
\item Riemann初始条件，$\rho$的计算结果等高线图像类似图\ref{rur}
\begin{figure}[h!]\centering
\includegraphics[width=14cm]{../Riemann-Laxwendrof-rou.eps}
\caption{Riemann初始条件的Lax-Wendrof格式计算解$\rho$}\label{rur}
\end{figure}
\end{enumerate}

收敛阶分析，对x轴\[-1,1\]分为40，80，160，320个网格计算了误差，作出下图。
\begin{enumerate}
\item $L^1$收敛阶，如图\ref{l1u}，双对数坐标下。
\begin{figure}[h!]\centering
\includegraphics[width=14cm]{../L1-convegence-order-of-lax-wendrof.eps}
\caption{网格大小与$L^1$误差的双对数图}\label{l1u}
\end{figure}
\item $L^2$收敛阶，如图\ref{l2u},双对数坐标。
\begin{figure}[h!]\centering
\includegraphics[width=14cm]{../L2-convegence-order-of-lax-wendrof.eps}
\caption{网格大小与$L^2$误差的双对数图}\label{l2u}
\end{figure}
\item $L^{\infty}$收敛阶，如图\ref{liu},双对数坐标。
\begin{figure}[h!]\centering
\includegraphics[width=14cm]{../L-infinity-convegence-order-of-lax-wendrof.eps}
\caption{网格大小与$L^{\infty}$误差的双对数图}\label{liu}
\end{figure}

\end{enumerate}
由图可知，Lax-wendrof格式是二阶收敛。

\begin{thebibliography}{99}
% adaptive
\bibitem{thz01}
李治平，《偏微分方程数值解讲义》，2009
\bibitem{thz02}
汤华中，《计算流体力学讲义》，2011
\end{thebibliography}

\end{document}
